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We analyze the methodology and the performance of subsystem density functional theory (DFT) 
with meta-generalized gradient approximation (meta-GGA) exchange-correlation functionals for 
non-bonded systems. Meta-GGA functionals depend on the Kohn-Sham kinetic energy density 
(KED), which is not known as an explicit functional of the density. Therefore, they cannot be 
directly applied in subsystem DFT calculations. We propose a Laplacian-level approximation to 
the KED which overcomes the problem and provides a simple and accurate way to apply meta- 
GGA exchange-correlation functionals in subsystem DFT calculations. The so obtained density and 
energy errors, with respect to the corresponding supermolecular calculations, are comparable with 
conventional approaches, depending almost exclusively on the approximations in the non-additive 
kinetic embedding term. An embedding energy error decomposition explains the accuracy of our 
method. 


I. INTRODUCTION 

Subsystem density functional theory is nowadays 
attracting increasing interest in the density functional 
theory (DFT) [3,@ community. This is due to its promise 
of achieving potentially exact results at a reduced com¬ 
putational cost, as well as to the high insight into the 
nature of interacting systems provided by the associated 
embedding potentials. Thus, numerous applications re¬ 
lated to non-covalent as well as covalent bonded 

systems [28l - [3lj | have been considered. In addition, the 
frozen density embedding (FDE) method n, m, m has 
emerged as a practical tool for efficient simulations of 
different properties [l^ . . We also recall that the 

FDE method with the iterative procedure of Ref. [3 is a 
computational implementation of the subsystem DFT. 

However, the accuracy of subsystem DFT calculations 
is practically limited by two factors. First, the term de¬ 
scribing the interaction energy between different subsys¬ 
tems depends on the non-additive kinetic energy (KE), 
which must be described by an explicit density func¬ 
tional. Second, the embedding potential, which is re¬ 
quired to describe the mutual interaction between the 
subsystems, must be a local multiplicative potential. 
Thus, it can contain only local or semilocal approxima¬ 
tions for the non-additive contribution of the kinetic and 
exchange-correlation (XC) terms to the embedding po¬ 
tential. Nevertheless, both limiting factors have currently 
been, at least partially, overcome. In fact, past years 
have seen the development of numerous KE functionals 
which can be suitable for subsystem DFT: GGA function¬ 
als (l^.[^. l39l - l4^ . Laplacian-level meta-GGA functionals 
and non-decomposable approac h For a recent 
review of all KE functionals see Ref. |^. Moreover, sev¬ 
eral works have extended the subsystem formulation of 


DFT beyond the conventional Kohn-Sham (KS) frame¬ 
work, considering e.g., hybrid functionals [d^, embedded 
interacting wave functions l45|. orbital-dependent effec¬ 
tive exact exchange methods [2lL , or density matrix 
[ 4 ^ . Nevertheless, to date, no attempt has been made 
to consider subsystem DFT calculations using meta gen¬ 
eralized gradient approximation (meta-GGA) XC func¬ 
tionals. 

A meta-GGA XC functional is defined by the general 
formula 

Exc = J e 3 ;c(p(r), Vp(r),r^®(r))dr , (1) 

where Cxc is the XC energy density, p is the electron 
density and 

-| occ. 

-'‘4r) = Ii:iWP(r)|" (2) 

is the positive-defined KS kinetic-energy density (KED), 
with being the occupied KS orbitals of the system. 
Meta-GGAs are attracting increasing popularity p§ - [^ 
because they can satisfy numerous exact constraints of 
the XC energy [a EiJsTIJpJp , achieve a remark¬ 
able level of accuracy [ , and describe 

excitonic effects in crystals |75|. In short, meta-GGA 
functionals has much larger accuracy/computational cost 
than GGAs, and thus should be preferred to the latter. 

However, is not an explicit functional of the den¬ 
sity, thus the implementation of meta-GGA functionals 
within the conventional KS scheme is not straightforward 
[ 73 , since it requires the calculation of /Sp. Thus, 
meta-GGAs are often implemented within a Generalized 
Kohn-Sham scheme (GKS) [t^- For this reason, so far, 
meta-GGA XC functionals have never been employed in 
subsystem DFT calculations. 
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In this paper we consider this issue and develop the 
theory and the methodology required to perform subsys¬ 
tem DFT calculations at the meta-GGA level. In par¬ 
ticular, we consider proper semilocal approximations for 
the XG embedding contributions and test them on a set 
of non-covalent complexes assessing the accuracy of the 
resulting energies and densities. 

Thus, the paper is organized as follows: in section [II] 
we present the general theory for subsystem DFT with 
arbitrary orbital-dependent XG functionals (thus includ¬ 
ing both meta-GGA and hybrid functionals) as well as 
different schemes to approximate the KED; Computa¬ 
tional details are reported in section Hill Results for non- 
covalent complexes are presented in section HVl Finally, 
in Section |V] "we summarize our conclusions. 


II. THEORY 


A. Subsystem DFT with arbitrary 
orbital-dependent XC functionals 


Within the subsystem formulation of density func¬ 
tional theory a given system is partitioned into two sub¬ 
systems A and B, each defined by its nuclear poten¬ 
tial fnuc(i’) respectively. Accordingly, the 

electron density of the total system is constructed as 
PA+B = PA + Pb, where the two subsystem densities 
integrate to Na and Nb, respectively. For simplicity we 
focus here on the case where Na and Nb are integer 
numbers. 

The ground-state solution of the problem is given by 
the set of coupled equations 


5Fbk[pa] 
SpA{r) 
^Fhk [pb] 
<5ps(r) 


= M (3) 

= M, (4) 


where Fhk is the universal functional of Hohenberg and 
Kohn, p, is the chemical potential (which is, in the case of 
the exact theory, equal in the two subsystems and equal 
to the supramolecular one [ 1 , [ISllll) 




SFbk [pi + Pii] 

Spi{r) 


5Fhk[pi] 

Spiir) 


-I- V 


II 

ext 


(r) 


(5) 


is the embedding potential with I = A, B and 
II = B,A, respectively (this convention will be used 
throughout). 

In this work we consider for Fhk the following parti¬ 
tion: 


FbkIp] = Fs[p] + J[p] (6) 

with 

(7) 


where p is any A^-representable electron density, J is the 
classical Goulomb energy, T is the KE operator, is a 
proper orbital-dependent XC functional (e.g. in the form 
given by Eq. ([T])), $ is a Slater determinant, and {cjii} are 
its single-particle orbitals. Equation © is quite general 
and includes not only meta-GGA functionals, but all the 
orbital-dependent ones. Then, Eqs. m and dH) become 


SFs[pi] 

Spi{r) 


6J[pi] 


+ = P , 


( 8 ) 


and the embedding potential of Eq. m can be written 
as 


^^emb(l’) = ^^ext(r) + 


SJjpii] 

hii{^) 


5Ff-^<^[pi,Pu] 

Spi{r) 


(9) 


where 


pb] = F4pa + pb] - Fs[pa] - F4pb] (10) 

is the non-additive contribution to the kinetic plus XC 
energy and we used the fact that the Coulomb potential 
is additive. 

At this point, following the GKS scheme Hi, [73 , we 
introduce, for each subsystem (for example I) an auxil¬ 
iary system of particles having the following properties: 
(i) it has the same ground-state density as our original 
embedded subsystem /; (ii) it is described by a single 
Slater determinant (iii) the ground-state energy is 
the minimum of the energy functional 


E[<I>^] = {<I>^\f\<I>^)+E,4{c^l}] + J w^{r)pi{r)d^r , 

( 11 ) 

where m is a (yet unknown) external local potential. 

The ground-state energy of this system is defined via 
the constrained search procedure as 


Eo 


min 

pi^Ni 


min 

pi^Ni 


|Cs[p/]+ J w^(r)p/(r)d3r| . 


( 12 ) 


Hence, the ground-state is described by the Euler equa¬ 
tion 


SFsjpi] 

Spi{r) 


+ w^(r) = p . 


(13) 


Comparing Eqs. ([SI) and da and making use of the 
property (i) we thus find 


=vLci'^) + 


^J]pi] 

Spi{r) 


^^Lb(r) ■ 


(14) 


On the other hand, properties (ii) and (iii) imply di¬ 
rectly that the ground-state of the auxiliary system is 
described by the set of single-particle equations 


--V^ + w\r) 




F,[p] = min {($|f|$) 


ef<('[(!•) , (15) 
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where we assumed real orbitals and ej are Lagrange mul¬ 
tipliers to ensure orbital orthonormality. For a meta- 
GGA functional the last term on the left hand side can 
be evaluated as 


SE, 


S(j>i 




(16) 

Combining Eq. (1141) with Eq. (fTSl) the operational equa¬ 
tions to solve the subsystem ground-state problem are, 
finally: 


2 


2 


A 

SpAir) 

1 SE,,[{^t}] 

2 <5<^f(r) 

SJ[p^ B 
SpB{r) ^ 


(r) 


.(r) 


= eUfir) 


(r) 


(r) 






(17) 


(ii) The functional depends explicitly on the or¬ 
bitals and has only an implicit dependence on the density. 
Therefore, to make the required functional derivatives 
special techniques are required such as the optimized 
effective potential method [90l - [^ . 

Consequently, in search of a practical computational 
procedure to perform subsystem DFT calculations with 
meta-GGA XC functionals, we propose to consider, in 
analogy with Refs. [U 13 in Eq. ra the semilocal 
approximation 

Fs[p] « Fs[p]=fs[p]+E^c[p], (20) 

_ F^‘^dd^p] = + E2f‘^[p] , ( 21 ) 

where the tilde denotes that an approximated (semilocal) 
functional of the density is used. 

For the kinetic term standard semilocal approxima¬ 
tions can be employed [IlillllSlIil. For the XC 
part, instead, two main possibilities can be envisaged: 


2 <5<^f(r) 


f</>f(r) 


(18) 


with the embedding potential given by Eq. (O. 

Note that, if the embedding potential is treated exactly, 
Eqs. (fT71) and (fTSl) admit in general multiple solutions 
d) 13 l82l|. Once the orbitals have been obtained, the 
total electron density is computed as 


Na Nb 

PA+B(r) = ^ |<())^(r)|%^ |(;if (r)l^ , (19) 

2 — 1 2—1 

whereas the total electronic energy is calculated as 

Ea+b[pa,Pb] = ($^|Ts|$^) + 

+-E'xc[{</>f}] + Exc[{4>f}] + J[pA + Pb] + 

+ J {PA{r) + PB{r)) (Vnuci^) + ^'fuc(l■)) pb] 


i) As a first simple option it is possible to use for E^c 
the GCA functional “most similar” to the meta- 
GGA functional used for s ubsy stems calculations. 
For example, if the TPSS [3l functional is used 
for subsystems calculations, the natural choice will 
be to use the PBE functional for the non-additive 
contribution. In fact, the TPSS functional has been 
constructed as an extension of the PBE functional. 
This choice resembles that of Refs. [2l|)ll3 and may 
be expected to yield reasonable results. The accu¬ 
racy of this combination will be verified in section 

El 

ii) A second, possibly better, choice is to retain for 

Exc the meta-GGA form, but replacing r with a 
suitable semilocal approximation. We recall that 
this is in general a very hard task However, 

for our purposes we only need the non-additive XC 
energy and, as shown in section m a large error 
cancellation effect can thus be expected. 


Note that the formalism introduced above is quite gen¬ 
eral and can be applied to GGA, hybrid-GGA, and meta- 
GGA functionals. 


C. Semilocal models for the kinetic energy density 


B. Non-additive embedding contributions 

To perform embedding calculations according to the 
theory detailed in Sec. Ill Al we need to compute the non¬ 
additive contribution (see Eq. (I^ l and 

its derivatives with respect to pA and ps (see Eq. ®)- 
This is not an easy task for two main reasons [dl] - 
(i) The computation of Fs[pA+ Pb] requires the knowl¬ 
edge of the ground-state Slater determinant of the total 
system, which is not available by definition in a 

subsystem calculation. It can only be obtained by an 
inverse KS [3 - [3| calculation starting from the ground- 
state total density pa a- b (for some examples within FDE 
see Refs. HH) 134^ . 


In this subsection we consider the construction of 
semilocal models for the KED. However, since the KED 
is not an observable, it is defined only up to a gauge in¬ 
tegrating to zero (and vanishing in the functional deriva¬ 
tive). Thus, to fix our working definition of the KED we 
decide to consider here only the positive-defined KED 
(Eq. (HI; see Ref. |93 for a discussion on this topic), 
which is also the most commonly used in meta-GGA XC 
functionals. 

To model the positive-defined KED in our subsystem 
DFT calculations we consider the following two semilocal 
approximations: 


T 

f 


TF 


1 

■ T = T 


20 


-r q = T 


( 22 ) 

(23) 
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where = ( 3 / 10 )( 37 r^)^/^p®'^^ is the Thomas- 

Fermi (TF) kinetic energy density mm, = 
(5/3)t"'’^s^ = |VjOp/( 8 p) is the von Weizsacker kinetic 
energy density |lOl| (with s = |Vp|/[ 2 ( 37 r^)^/^p'^/^] the 
reduced gradient), 7 -revAPBEK jg revAPBEK kinetic 
energy density [3 [13; and q = V^/ 9 /[ 4 ( 37 r^)^/^p^/^] is 
the reduced Laplacian. 

The model of Eq. is a simple GGA model 

that is exact for a uniform density perturbed by a small- 
amplitude short-wavelength density wave and motivated 
by the basic requirements that r « in the slowly- 
varying density limit and r Ki |l 02 j| in tail regions 
(where 0) and iso-orbital regions. Moreover, this 

simple model fulfills the important constraint < r, 
i.e. that 0 < z < 1, with z = r'^/r, which has actually 
been used in the construction of several meta-GGA XG 
functionals [48l |. 

The model of Eq. (1251) is a Laplacian-level meta- 
GGA model and requires several considerations: 

i) We first recall that (20/9)r'^^q = (l/ 6 )V^p. Thus, 
the kinetic energy and potential corresponding to 

are identical with the revAPBEK ones (recall 
that a term proportional to the Laplacian of the 
density does not contribute to the energy and the 
potential). If the revAPBEK functional is used for 
f, and in Eqs. ([201) and (EB), then the re¬ 

vAPBEK KE approximation is “de facto” the only 
functional approximation used in the subsystem 
DFT meta-GGA calculation. In fact, in our imple¬ 
mentation we use the same KE functional in both 
and 

ii) The term containing q is of fundamental impor- 

tance to reproduce the correct KE density [^[^). 
In Eq. (1251) the coefficient of the reduced Laplacian 
term comes from the lowest-order Laplacian con- 
tributio n to the s econd-order gradient expansion of 
the KE [l0l[^. Other coefficients could be used 
as well but we found that the non-empirical 

one in Eq. (1251) is quite accurate for our purpose, 
even if Eq. (1231) is not exact in the asymptotic re¬ 
gion (see Appendix lAI) . 

iii) The revAPBEK functional, which is used in the 
definition of r^, recovers by construction the 
mo dified second-order gradient expansion of the 
KE [l05j . which was constructed from semiclassical 
theory of atoms. Actually, the revAPBEK func¬ 
tional, which does not contain any empirical pa¬ 
rameter fitted on kinetic energies, has been found to 
be very accurate for the description of non-covalent 
complexes within subsystem DFT [ 23 . However, 
the accuracy of the model is not mandatorily 
related to the use of revAPBEK: other state-of- 
the-art GGA [H, HO, Ili-Iil and meta-GGA [H] 
KE approximations can be expected to yield simi¬ 
lar accuracy. 



distance [bohr] 


FIG. 1: Plot of z = r'^/r (panel a), a = (r — 

(panel b) and (panel c) for the Ne dimer along the 

main axis computed using the kinetic energy density models 
of Eqs. (I 22 I) and (l23l) . 

iv) The Laplacian term will diverge at the core {q —> 
— 00 ). However, such a bad behavior in the core 
region is not a problem for subsystem DFT calcu¬ 
lations, since for the non-additive XC energy and 
potential these contributions cancel almost com¬ 
pletely This cancellation is shown in Figure 6 of 
Ref. [2^ where a reduced gradient and Laplacian 
decomposition of the non-additive KE is reported, 
showing that only small values of q contribute sig¬ 
nificantly. 

To show the physical significance of the two models 
introduced above and offer a preliminary test of the ex¬ 
pectable performance, we consider their application to 
the test case of the Ne dimer, an example for weakly in¬ 
teracting systems. As meta-GGA exchange correlation 
functional we consider the TPSS [13 one. This is in fact 
one of the first and most popular non-empirical meta- 
GGAs and is used here as an exemplary case. 

In Fig. [TJi we report the quantity 2 = r^/r for the 
two KE density approximations Eqs. (1221) and (1251) as a 
function of the distance d along the dimer axis. In the 
middle of the bond {d = 0) all curves go to zero, because 
the gradient of the density and thus vanish. The 
model is accurate only at few points: |d| si I bohr and 



















5 


near the core. However, exa ctly a t the nucleus position 
the exact /t is close to 1 [l02j| : on the other hand at 
this point the density is very large so that ^ , 

thus jT^ is much less than 1. The Laplacian model 
is everywhere (but in the core) much more accurate 
than r^. In particular, it reproduces almost exactly 
in the region \d\ < 1.5 bohr. 

Alternatively, in Fig. [Hr we report the quantity 



which is the main, direct ingredient of the meta-GGA 
(TPSS) exchange enhancement factor and thus on the 
calculation of the exchange energy. Differences about 
the various approaches are now more evident. For the 
model 0 = 1 everywhere in the space (by construction): 
this can be a good approximation only for slowly vary¬ 
ing densities. Instead, for molecular systems, the exact 
alpha shows significant oscillations and it is very high 
(a ss 17.6) in the bond. These oscillations are correctly 
reproduced by the model, which yields a very accu¬ 
rate value at the bond (a « 16.0). is significantly 
different from the exact one only near the core, where it 
is actually negative. Finally, In Fig. [TJ: we report the 
quantity 

AFJP®® = FjP®®(p, Vp, f) - FjP®®(p, Vp, r^®) (25) 

where is the TPSS exchange enhancement factor 

[ 4 ^ . so that 

where is the LDA exchange energy per particle. The 
quantity indicates how much the approximation 

in T will impact on the accuracy of the exchange energy. 

The plots in Fig. [T]3 confirm the high accuracy of the 
model, whereas the model shows quite larger dif¬ 
ferences in the region |d| < 0.5 au. 

In section 113 we will consider the performance of the 
two models for subsystem DFT calculations. 

III. COMPUTATIONAL DETAILS 

To assess the possibility of performing subsystem DFT 
calculations using meta-GGA functionals we carried on 
test simulations on different non-covalent complexes. To 
this end, for simplicity, we considered as meta-GGA XG 
functional the TPSS one [i^. Other meta-GGA XG func¬ 
tionals will be considered in detail in a future publication. 

In our calculations we used different approximations 
for the non-additive XC terms. We refer to each of these 
using the notation methodI/method2, that denotes that 
methodl was used for the subsystems and method2 to 
compute the non-additive XG contribution. In more de¬ 
tails: 


1. As a simple choice we computed the non-additive 

XC contributions using the PBE XC functional 
|l06j| . i.e. we set E^cc = . This approach 

is labeled TPSS/PBE. 

2. Alternatively, we computed the non-additive XC 
terms using the XC TPSS functional but using the 

model of Eq. to mimic the positive-defined 
KED. This approach is named TPSS/TPSS-I. 

3. Finally, we considered the same case as above but 
using the model of Eq. (l23l) . This approach is 
labeled TPSS/TPSS-L. Note that in this case, be¬ 
cause of the negative divergence of q, the model we 
use for r is not guaranteed to respect the bound 

< T. Nevertheless, the TPSS functional is nu¬ 
merically well defined also for z < 0 or z > 1; more¬ 
over these values occur only near the core which 
is not relevant for non-additive contributions, see 
point iv) of section HTCI 

We remark that the first two approximations are GGA 
ones, while the last one is a Laplacian-level meta-GGA 
method. For compa rison also subsy stem DFT calcula¬ 
tions using the PBE |l06j | and PBEO [l07l Il08j| XC func¬ 
tionals were considered. The former one requires no ap¬ 
proximations for the non-additive XC terms and imple¬ 
ments the PBE/PBE approach; the latter one instead 
uses a semilocal approximation as described in Ref. |3 
and yields the PBEO/PBE approach. 

In all calculations the non-additive kinetic contribu¬ 
tions were computed using the revAPBEK kinetic func- 
tiona]_H9jM and a supermolecular def2-TZVPPD basis 
set [I09l IllOjl was employed. As the aim of this work is 
not to verify the absolute accuracy of the embedding ap¬ 
proach (which depends critically on the KE approxima¬ 
tion) , but if the additional errors due to the non-additive 
XC approximation can be reduced or not, we believe that 
checking one (accurate) KE functional is enough. 

All calculations have been performed usin g th e FDE 
script [15 of the TURBOMOLE program package |lll| . The 
calculation of the matrix elements of the non-additive 
TPSS-L XC functional (which is a Laplacian-level meta- 
GGA functional), has been performed as described in the 
Appendix A of Ref. . 

The complexes considered for the tests have been di¬ 
vided into four groups according to the character domi¬ 
nating their interaction: 

- WI (weak interaction): He-Ne, He-Ar, Ne-Ne, Ne- 
Ar, CH 4 -Ne, CgHe-Ne, CH 4 -CH 4 ; 

- DI (dipole-dipole interaction): H 2 S-H 2 S, HCl-HCl, 
H 2 S-HCI, CH 3 CI-HCI, CH 3 SH-HCN, CH 3 SH-HCI; 

- HB (hydrogen bond): NH 3 -NH 3 , HF-HF, H 2 O- 
H 2 O, HF-HCN, (HC0NH2)2, (HC00H)2; 

- CT (charge transfer): NF 3 -HCN,C 2 H 4 -F 2 ,NF 3 - 
HCN, C 2 H 4 -CI 2 , NH 3 -F 2 , NH 3 -CIF, NF 3 -HF, 
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C 2 H 2 -CIF, HCN-CIF, NH 3 -CI 2 , H 2 O-CIF, NH 3 - 
CIF. 

The reference geometries a nd binding energies were taken 
from Refs, i, |M Ell, IHl 

The error on the total embedding energy was computed 
as the difference between the energy obtained from a sub¬ 
system DFT calculation (i.e. Eq. (l20l) l and the energy 
(^GKS) obtained from the corresponding supermolecular 
conventional calculation 1 ^ , i.e. 

= Ea+b[pa, Pb] - E^^^iPGKs] (27) 

where pA and pB are approximated embedded subsys¬ 
tems densities, due to the approximations in Eq. (EU). 

The performance of the different approaches was eval¬ 
uated, within each group of complexes, by computing the 
mean absolute error (MAE) and the mean absolute rel¬ 
ative error (MARE) with respect to reference binding 
energies [l^. Instead, to assess the performance of the 
methods for all the different classes of systems, we con¬ 
sidered the quantities 



FIG. 2: Plot of - Appbe/pbe{z) (for the Ne 

dimer along the main axis) computed using the kinetic energy 
density models of Eqs. ([^ and (1^ as well as the PBE XC 
functional for the non-additive XC embedding potential. 


IV. RESULTS 


rwMAE 


rwMARE 


1 

4 

1 

4 


E 

i=WI,DI,HB.CT 


E 

i=WI,DI,HB,CT 


/ MAE, \ 
\{MAE,)) 

/ MARE, \ 
[{MARE,}) 


(28) 

(29) 


where {MAEi) {{MAREi)) is the average MAE (MARE) 
among the different methods considered for the class of 
systems i. In this way, all the different classes of systems 
will have the same influence of the rwMAE (rwMARE) 
and methods with rwMAE (rwMARE) smaller than I 
will have better performance than the average. 

The errors on the embedding densities were studied by 
considering the deformation density 

Ap{r) = pA{r) + pB{r) - , (30) 

where denotes the density obtained from a con¬ 

ventional GKS calculation. A quantitative measurement 
of the absolute error associated with a given embedding 
density was then obtained by computing the embedding 
density error 


\Ap{r)\dr, (31) 

with N the number of electrons. In the evaluation of 
f only valence electron densities were considered. Core 
densities are in fact much higher than valence ones and 
would largely dominate. On the other hand, core den¬ 
sities are not very important for the determination of 
chemical and physical properties of the interaction be¬ 
tween the subsystems, which are of interest here. The 
performance of the different approaches was evaluated 
by computing the MAE and the rwMAE. 


The errors f on the embedding densities for the dif¬ 
ferent methods are reported in Table HI The best per¬ 
formance is observed for the PBEO/PBE method, which 
gives the smallest error for all the systems investigated. 
As explained in Refs. [2l|, [U, this fact traces back 
to the reduced self-interaction error of the PBEO func¬ 
tional, which reduces the overlap between the subsys¬ 
tem densities. All other methods yield very close re¬ 
sults with a rwMAE in the range 0.97 — 1.14. Actually 
the TPSS/TPSS-L has the lowest rwMAE between these 
ones, showing that meta-GGA subsystem calculations 
can perform even better than conventional GGA calcula¬ 
tions, despite the former include an additional approxi¬ 
mation. Among meta-GGA methods, the TPSS/TPSS-L 
approach is the best for WI systems (MAE=0.13) and for 
DI (MAE=2.46), whereas TPSS/PBE is the best for HB 
and GT. 

The integrated measure f however provides only an ab¬ 
solute measure of the error on embedding densities, but 
cannot tell anything on how the error on the density is 
distributed in the space and what are the roles of the ki¬ 
netic and XG approximations to determine such an error. 
Here, we aim at understanding better the importance 
of different approximations used in the non-additive XG 
term of the embedding potential. Thus, we consider in 
Fig. [^the plot along the bond axis of the Ne-Ne complex 
(taken as an example) of the plane-averaged XG defor¬ 
mation density Ap'^‘^*^°‘^{z) —(z), where the 
plane-averaged deformation density is 

Ap{z) = j \Ap{r)\dxdy . (32) 

This quantity provides in fact a measure, point by point, 
of the embedding density error due to the XC approxima¬ 
tion: the PBE/PBE is in fact taken as reference because 
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TABLE I: Errors on the embedding densities for different methods as obtained using Eq. (ED). At the bottom of each group 
of results the mean absolute error (MAE) is reported. The last row report the global rwMAE (see text for details). The best 
result of each line is highlighted in bold style. A star indicates the best result among meta-GGA methods. 


Complex 

PBE/PBE 

PBEO/PBE 

TPSS/PBE 

TPSS/TPSS-1 

TPSS/TPSS-L 



weak interaction (WI) 



He-Ne 

0.05 

0.02 

0.05* 

0.05* 

0.05* 

He-Ar 

0.06 

0.06 

0.12 

0.08 

0.05* 

Ne-Ne 

0.04 

0.02 

0.06 

0.04 

0.03* 

Ne-Ar 

0.06 

0.04 

0.13 

0.08 

0.05* 

CH 4 -Ne 

0.07 

0.05 

0.16 

0.13 

0.06* 

CgHe-Ne 

0.13 

0.11 

0.25 

0.27 

0.13* 

CH 4 -CH 4 

0.60 

0.50 

0.84 

0.79 

0.53* 

MAE 

0.14 

0.11 

0.23 

0.21 

0.13* 



dipole-dipole interaction (DI) 



H 2 S-H 2 S 

1.85 

1.62 

1.89 

1.81 

1.70* 

HCl-HCl 

1.87 

1.49 

1.88 

1.91 

1.75 

H 2 S-HCI 

3.70 

2.97 

3.59 

3.74 

3.56 

CH 3 CI-HCI 

2.38 

1.91 

2.36 

2.40 

2.24 

CH 3 SH-HCN 

1.72 

1.61 

1.74 

1.64 

1.58 

CH 3 SH-HCI 

4.08 

3.32 

3.90 

4.12 

3.95 

MAE 

2.60 

2.15 

2.56 

2.60 

2.46* 



hydrogen bond (HB) 



NH 3 -NH 3 

1.79 

1.67 

1.87 

1.85 

1.74 

HF-HF 

1.53 

1.19 

1.56 

1.62 

1.50* 

H 2 O-H 2 O 

2.01 

1.72 

2.05 

2.11 

* 

00 

NH 3 -H 2 O 

3.11 

2.69 

3.06* 

3.19 

3.08 

HF-HCN 

2.77 

2.38 

2.62* 

2.84 

2.75 

(HC 0 NH 2)2 

2.72 

2.49 

2.71* 

2.76 

2.65 

(HC00H)2 

3.35 

2.94 

3.23* 

3.45 

3.37 

MAE 

2.47 

2.15 

2.44* 

2.55 

2.54 



charge transfer (CT) 



NF 3 -HCN 

0.29 

0.24 

0.40 

0.43 

0.26* 

C 2 H 4 -F 2 

6.35 

2.75 

5.68* 

5.79 

5.79 

NF 3 -HNC 

0.58 

0.49 

0.58 

0.63 

0.55* 

C 2 H 4 -CI 2 

5.77 

4.32 

5.85* 

6.08 

6.28 

NH 3 -F 2 

9.60 

4.38 

8.48 

8.60 

8.58* 

NF 3 -CIF 

1.73 

0.99 

1.59 

1.54* 

1.68 

NF 3 -HF 

0.95 

0.75 

0.91 

0 . 88 * 

0.91 

C 2 H 2 -CIF 

6.02 

4.32 

5.97* 

6.77 

6.51 

HCN-CIF 

3.21 

2.33 

3.08* 

3.40 

3.30 

NH 3 -CI 2 

7.60 

5.48 

7.42* 

8.25 

8.06 

H 2 O-CIF 

5.06 

3.42 

4.98 

5.54 

5.39* 

NH 3 -CIF 

11.19 

9.37 

11 . 00 * 

12.37 

12.06 

MAE 

4.86 

3.24 

4.66* 

5.02 

4.95 

rwMAE 

1.00 

0.79 

1.12 

1.14 

0.97* 








it includes only approximation due to KE. We see that, 
in agreement with the results of Fig. [H the TPSS-L ap¬ 
proximation performs very well, introducing only very 
small errors in the calculation of the embedding density. 
A larger effect is instead obtained for TPSS-1, while the 
use of the PBE XC functional to compute the embedding 
potential yields considerably larger differences. 

We now turn to discuss the embedding energy errors, 
which are reported in Tab. m In this case the hy¬ 
brid PBEO/PBE is not the best for all systems, as it 
was in Tab. H] In fact, for dipole-dipole and hydrogen 
bond complexes the TPSS/TPSS-L method is the most 
accurate, closely followed by the PBE/PBE approach, 
whereas the PBEO/PBE is not so accurate for these sys¬ 
tems [m [ 13 . On the other hand, PBEO/PBE is very 
accurate for weakly-interacting and charge-transfer sys¬ 
tems 

A more detailed discussion of the trends is provided 
in subsection IIV Al where we perform an energy de¬ 
composition analysis of the embedding energy errors. 
Here we just note that, concerning the meta-GGA ap¬ 
proaches, our simple Laplacian-level meta-GGA approx¬ 
imation (TPSS/TPSS-L) is significantly more accurate 
than the GGA TPSS/TPSS-1 and TPSS/PBE ones and 
can also slightly outperform the use of a simple GGA 
XG functional, such as PBE. In fact, for TPSS/TPSS-L, 
PBE/PBE, TPSS/PBE, and TPSS/TPSS-I we find over¬ 
all rwMAEs of 0.84, 0.93 and 1.21, 1.24, respectively. 


A. Error decomposition analysis 

The embedding energy error for meta-GGA as well as 
for hybrid functionals depends on two distinct approx¬ 
imations, the KE and the XC energy. As discussed in 
Ref. [12 this causes a subtle error cancellation effect. To 
understand better the error compensation issue in meta- 
GGA subsystem DFT calculations and analyze in detail 
the sources of different errors, we perform in the following 
an error decomposition analysis. 

Following Ref. [12 we thus write the error on the em¬ 
bedding energy as 


with 

AE — AT, -|- AD AExc j 

(33) 

at. 


(34) 


- (r.[{4^l}] - TslW^}] - T,[{</^}]^ 

) 

AD 

= EextiPA+s] + J[PA+b] + Exc[PA+b] ~ 
-{Ecxt[p^^^]+J[p^^^]+Exc[p^^^]) 

(35) 

AExc 

= Exc[p^^^] - Exc[pa] - ExcIpb] - 

(36) 


- {ExciWZl}] - ExclW^}] - ExciW^}]) . 


Here the first term (AT^) describes the error associated 
with the KE approximation, the second term (AH) is a 


relaxation term directly related to the embedding density 
error, and the third one (AE^c) measures the error due to 
the XC approximation. This latter term will be negative 
when the approximate XC functional overestimates with 
respect to the non-local one and positive in the opposite 
case. In analogy to the previous notation, (pA and (j)B 
denote single particle KS orbitals of subsystems A and B, 
respectively, as obtained from approximated embedding 
calculations. 

The results of the energy decomposition analysis are 
reported in Table m The contributions due to AT^ and 
AD are reported summed together because both terms 
yield very large values but with opposite sign, thus they 
only contribute to the total error through a strong error 
cancellation. 

For each group of complexes as well as for the over¬ 
all test set we report, for each component of the energy 
decomposition, the MAE and the MARE. Moreover, for 
AExc we report the XC differential error (XCDE) and 
the XC differential relative error (XCDRE), defined as 


XCDE 

XCDRE 


1 

-^|AH,|-|AT,,, + AA| (37) 

1 ^|AE,|-|AT,,,, + AD,l 


where the sums extend over all the N systems in the set. 
A positive value of these statistical indicators denotes 
that the XC approximation has a bad effect on the abso¬ 
lute total embedding energy error, increasing it. On the 
contrary, a negative value indicates that the employed 
XC approximation reduces the error (presumably by er¬ 
ror cancellation). 

Inspection of the table shows that the ATg + AD val¬ 
ues are very similar for all the methods. This reflects 
the fact that we used the same KE approximation in 
all calculations and that in all cases the final embed¬ 
ding densities are quite similar. On the other hand, the 
values of AExc are more differentiated between the var¬ 
ious methods. In particular, much smaller values are 
found in general for the TPSS/TPSS-L method (global 
MARE 4%). The TPSS/TPSS-1 (global MARE 25%) 
and TPSS/PBE (global MARE 32%) accuracy is much 
less. This confirms that the TPSS/TPSS-L benefits of a 
much better XC approximation than the latter. 

Valuable information is also obtained by the inspec¬ 
tion of the XCDRE indicators. For WI and CT systems 
the TPSS/PBE and TPSS/TPSS-1 methods have neg¬ 
ative values of XCDRE. Thus the additional error due 
to the XC approximation yields (due to error cancella¬ 
tion) better total energies. This explains the results in 
Tab. n where TPSS/PBE (and TPSS/TPSS-1) shows 
a good accuracy for these systems. On the other hand, 
for DI and HB systems, the XCDRE values are positive, 
i.e. the additional error due to the XC approximation 
reduces the accuracy of the embedding energy. In these 
cases no error cancellation occurs and in fact TPSS/PBE 
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TABLE II: E mbeddin g: energy errors (mHa) for different methods and complexes. Accurate reference binding energies (Eb) 
from Refs. Ill2l . 1 1131 are also reported in the second column. At the bottom of each group of results the mean absolute error 


(MAE) and the mean absolute relative error 
star indicates the best result among the ones 

(MARE) are reported. The best result of each line is highlighted 
with the TPSS functional. 

in bold style. A 

Complex 

Eb 

PBE/PBE 

PBEO/PBE TPSS/PBE 

TPSS/TPSS-1 

TPSS/TPSS-L 

He-Ne 

0.06 

0.08 

weak interaction (WI) 

0.03 

0.03* 

0.06 

0.08 

He-Ar 

0.10 

0.05 

0.00 

- 0 . 01 * 

0.04 

0.06 

Ne-Ne 

0.13 

0.14 

0.06 

0.02* 

0.10 

0.13 

Ne-Ar 

0.21 

0.11 

0.04 

-0.04* 

0.06 

0.11 

CH 4 -Ne 

0.35 

0.12 

0.04 

-0.04* 

0.06 

0.12 

CgHe-Ne 

0.75 

-0.03 

- 0.10 

-0.51 

-0.25 

-0.01* 

CH 4 -CH 4 

0.81 

-0.38 

-0.41 

-0.82 

-0.54 

-0.27* 

MAE 


0.13 

0.10 

0.21 

0.16 

0 . 11 * 

MARE 


61% 

27% 

39%* 

52% 

59% 

H 2 S-H 2 S 

2.63 

-0.47 

dipole-dipole (DI) 
-0.84 

-1.16 

-1.07 

-0.49* 

HCl-HCl 

3.20 

0.07 

-0.37 

-0.70 

-0.62 

-0.02* 

H 2 S-HCI 

5.34 

0.40 

-0.42 

-0.54 

-0.71 

0.29* 

CH 3 CI-HCI 

5.66 

0.02 

-0.59 

-1.14 

-1.27 

-0.05* 

CH 3 SH-HCN 

5.72 

-1.18 

-1.57 

-2.09 

-2.04 

-1.02* 

CH 3 SH-HCI 

6.63 

0.73 

-0.34 

-0.64 

-1.06 

0.54* 

MAE 


0.48 

0.69 

1.05 

1.13 

0.40* 

MARE 


10 % 

16% 

24% 

25% 

9%* 

NH 3 -NH 3 

5.02 

-0.95 

hydrogen bond (HB) 
-1.32 

-1.69 

-1.63 

-0.80* 

HF-HF 

7.28 

0.79 

0.19 

-0.13 

-0.03* 

0.78 

H 2 O-H 2 O 

7.92 

- 0.20 

-0.79 

- 1.11 

-1.15 

-0.15* 

NH 3 -H 2 O 

10.21 

-0.44 

-1.28 

-1.47 

-1.75 

-0.36* 

HF-HCN 

11.33 

0.43 

-0.56 

-0.72 

-1.06 

0.49* 

(HC 0 NH 2)2 

23.81 

-4.21 

-5.30 

-5.95 

-6.87 

-3.42* 

(HC00H)2 

25.74 

- 1.88 

-3.69 

-3.94 

-5.61 

-1.37* 

MAE 


1.27 

1.88 

2.14 

2.59 

1.05* 

MARE 


9% 

13% 

16% 

18% 

8%* 

NF 3 -HCN 

1.67 

-0.41 

charge transfer (CT) 

-0.43 

-0.95 

- 0.88 

-0.31* 

C 2 H 4 -F 2 

1.69 

4.27 

1.92 

3.13* 

3.42 

3.87 

NF 3 -HNC 

2.31 

-0.13 

-0.51 

-0.78 

- 1.11 

-0.02* 

C 2 H 4 -CI 2 

2.60 

1.52 

-0.42 

0.30* 

-1.87 

1.70 

NH 3 -F 2 

2.88 

6.90 

2.98 

5.17* 

5.47 

6.07 

NF 3 -CIF 

2.92 

2.14 

0.82 

0.88 

0.15* 

1.95 

NF 3 -HF 

2.92 

0.91 

0.22 

0.05* 

-0.57 

0.86 

C 2 H 2 -CIF 

6.07 

3.71 

1.52 

2.40 

1.64* 

3.77 

HCN-CIF 

7.74 

1.62 

0.03 

0.28 

-0.27* 

1.51 

NH 3 -CI 2 

7.78 

2.84 

0.21 

1.64 

0.85* 

2.94 

H 2 O-CIF 

8.54 

2.42 

0.45 

1.17 

0.55* 

2.45 

NH 3 -CIF 

16.92 

4.44 

1.31 

2.35 

-0.33* 

5.75 

MAE 


2.61 

0.90 

1.59 

1.43* 

2.60 

MARE 


72% 

30% 

49%* 

53% 

67% 


rwMAE 

0.93 

0.79 

1.24 

1.21 

0.84* 

rwMARE 

0.98% 

0.78% 

1 . 10 % 

1.24% 

0.91%* 
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TABLE III: Embedding energy error decomposition (mHa) for different methods and complexes. At the bottom of each group 
of results the mean absolute error (MAE) and the mean absolute relative error (MARE) are reported. In addition for AE^c 
also the XC differential error (XCDE) and the XC differential relative error (XCDRE) are listed. 


Complex 

TPSS/PBE 

TPSS/TPSS-1 

TPSS/TPSS-L 

ATs + AD 

/\.Exc 

ATs + AD 

A£/xc 

ATs + ad 

^Exc 




weak interaction 




He-Ne 

0.08 

-0.05 

0.08 

- 0.02 

0.08 

- 0.01 

He-Ar 

0.05 

-0.06 

0.05 

- 0.02 

0.05 

0.00 

Ne-Ne 

0.13 

- 0.11 

0.13 

-0.04 

0.14 

- 0.01 

Ne-Ar 

0.10 

-0.14 

0.11 

-0.05 

0.12 

0.00 

CH 4 -Ne 

0.10 

-0.14 

0.11 

-0.05 

0.11 

0.00 

CeHg-Ne 

-0.06 

-0.46 

-0.05 

- 0.20 

- 0.01 

0.00 

CH 4 -CH 4 

-0.43 

-0.39 

-0.42 

- 0.12 

-0.38 

0.11 

MAE 

0.14 

0.19 

0.14 

0.07 

0.13 

0.02 

MARE 

60% 

63% 

61% 

23% 

61% 

5% 

XCDE/XCDRE 


-l-0.08/-21% 


-|- 0 . 02 /- 12 % 


-0.02/-5% 




dipole-dipole 




H 2 S-H 2 S 

-0.48 

- 0.68 

- 0.45 

-0.62 

-0.28 

- 0.21 

HCl-HCl 

0.06 

-0.76 

0.08 

-0.71 

0.19 

-0.16 

H 2 S-HCI 

0.38 

-0.93 

0.44 

-1.15 

0.67 

-0.38 

CH 3 CI-HCI 

0.01 

-1.15 

0.02 

-1.29 

0.22 

-0.28 

CH 3 SH-HCN 

- 1.20 

-0.90 

- 1.21 

-0.84 

-0.99 

-0.03 

CH 3 SH-HCI 

0.71 

-1.35 

0.73 

-1.80 

1.20 

-0.67 

MAE 

0.47 

0.96 

0.49 

1.07 

0.59 

0.29 

MARE 

10 % 

21 % 

10 % 

22 % 

11 % 

6 % 

XCDE/XCDRE 


-l-0.58/+14% 


+0.65/4-15% 


-0.19/-3% 




hydrogen bond 




NH 3 -NH 3 

-0.96 

-0.74 

-0.96 

-0.67 

- 0.88 

0.08 

HF-HF 

0.79 

-0.92 

0.79 

-0.82 

0.89 

- 0.11 

H 2 O-H 2 O 

- 0.21 

-0.90 

- 0.20 

-0.95 

- 0.02 

-0.13 

NH 3 -H 2 O 

-0.46 

- 1.01 

-0.45 

-1.30 

-0.13 

-0.23 

HF-HCN 

0.41 

-1.13 

0.44 

-1.51 

0.68 

-0.19 

(HC0NH2)2 

-4.20 

-1.74 

-4.36 

-2.52 

-3.49 

0.07 

(HC00H)2 

-1.87 

-2.07 

-2.04 

-3.57 

-1.31 

-0.06 

MAE 

1.27 

1.22 

1.32 

1.62 

1.06 

0.12 

MARE 

9% 

11 % 

10 % 

12 % 

8 % 

1 % 

XCDE/XCDRE 


4-0.87/4-6% 


+1.27/+8% 


+ 0 . 00 /+ 0 % 




charge transfer 




NF 3 -HCN 

-0.46 

-0.50 

-0.44 

-0.19 

-0.36 

0.04 

C 2 H 4 -F 2 

4.23 

- 1.10 

4.21 

-0.79 

3.83 

0.04 

NF 3 -HNC 

-0.14 

-0.63 

-0.13 

-0.49 

0.00 

- 0.02 

C 2 H 4 -CI 2 

1.51 

- 1.21 

1.55 

-2.05 

1.66 

0.04 

NH 3 -F 2 

6.84 

-1.67 

6.94 

-1.48 

6.38 

-0.31 

NF 3 -CIF 

2.11 

- 1.22 

2.17 

-1.25 

2.11 

-0.17 

NF 3 -HF 

0.89 

-0.84 

0.92 

-0.83 

1.03 

-0.16 

C 2 H 2 -CIF 

3.71 

-1.31 

3.84 

- 2.21 

3.79 

- 0.02 

HCN-CIF 

1.60 

-1.32 

1.68 

-1.95 

1.63 

- 0.12 

NH 3 -CI 2 

2.84 

-1.19 

2.90 

-2.05 

3.07 

-0.14 

H 2 O-CIF 

2.41 

-1.24 

2.50 

-1.94 

2.54 

-0.09 

NH 3 -CIF 

4.46 

- 2.11 

4.27 

-4.60 

5.88 

-0.13 

MAE 

2.60 

1.20 

2.63 

1.65 

2.69 

0.11 

MARE 

71% 

32% 

72% 

35% 

69% 

3% 

XCDE/XCDRE 


- 1 . 01 /- 22 % 


-1.36/-25% 


-0.09/-2% 

MAE 

MARE 

XCDE/XCDRE 

1.37 

44% 

0.94 

32% 

-0.06/-9% 

1.40 

44% 

1.19 

25% 

-0.11/-7% 

1.38 

43% 

0.13 

4% 

-0.07/-3% 
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TABLE IV: Components of the embedding energy error de¬ 
composition {AExc and ATg+AD; mHa), XCAR (see text for 
details), and embedding density errors (5), from TPSS/TPSS- 
L calculations (last column report in parenthesis the value of ^ 
for PBE/PBE) for various test complexes as a function of the 
intermolecular distance. Rq denotes the equilibrium distance 
in the complexes (Ne 2 , Ro = 3.09lA; HF-NCH, Ro = l.SOSA; 
(HC1)2, Ro = 3.872A). 


System 

R/Ro AU„c ATs + AD XCAR 



Ne2 

0.8 

0.02 

-0.17 

10 % 

0.34 

(0.34) 


1.0 

-0.01 

0.14 

6 % 

0.03 

(0.04) 


1.5 

0.00 

0.00 

0 % 

0.00 

(0.00) 

(HC1)2 

0.8 

-1.55 

-3.40 

31% 

9.31 

(9.82) 


1.0 

-0.16 

0.19 

45% 

1.75 

(1.87) 


1.5 

0.00 

0.03 

3% 

0.05 

(0.05) 

HF-NCH 0.8 

-0.58 

-1.81 

24% 

5.49 

(5.48) 


1.0 

-0.19 

0.68 

22 % 

2.75 

(2.77) 


1.5 

-0.02 

0.39 

5% 

0.32 

(0.37) 


and TPSS/TPSS-1 yield quite bad total energy (see Tab. 

CD). 

On the other hand, in the TPSS/TPSS-L method the 
XCDRE values are very small (and negative), showing 
the smallest error compensation in relation to the XC 
approximation. Note, in addition, that for all the con¬ 
sidered TPSS subsystem DPT calculations the AE^c 
values are comparable or smaller than for the hybrid 
PBEO/PBE method (see Table II of Ref. and note 
that these values include a prefactor 0.25). 

Finally, we note that the good accuracy of the 
TPSS/TPSS-L approach is maintained also for larger 
subsystems’ density overlaps. This is shown in Tab. IIVI 
where we report the AE^c, ATs + AD, and the density 
error computed with the TPSS/TPSS-L approach for 
several complexes at different intermolecular distances 
(smaller distances correspond to higher overlaps). Tab. 
CYlalso reports the XC absolute ratio (XCAR), defined as 
XCAR = \AExc\/ (I ATg-l-AD|-|- 1 Aifa;c|), which is an ab¬ 
solute (i.e. without error compensation) measure of the 
non-additive XC contribution to the total embedding en¬ 
ergy error. The data reported in Table IIVI clearly show 
that at shortest distance XCAR is not largely increased 
(as it happens instead for ^), but remains constant and 
in some cases it is even reduced. 


V. CONCLUSIONS 

Using the generalized Kohn-Sham framework, we ex¬ 
tended the subsystem DFT formalism to the use of meta- 
GCA functionals. For a practical application of the 
method we proposed several semilocal approximations for 
the non-additive XC energy. Two of these are based on 
simple models for the KED. 

The results of our subsystem DFT calculations show 


that all the proposed methods work reasonably well, 
giving density and energy embedding errors compara¬ 
ble with conventional calculations and close to hybrid 
subsystem DFT results. Nevertheless, a more detailed 
analysis shows that for the simplest approaches the fi¬ 
nal performance is the result of a quite significant error 
compensation. This effect is reduced only when more so¬ 
phisticated approximations for the non-additive XC term 
are used. In this respect we showed that this goal can 
be pursued by considering Laplacian-level meta-GCA 
approximations. Anyway, we remark that our TPSS- 
L approximation, despite giving promising results, is 
only a simple model used here to investigate the power 
of Laplacian-dependent approximations and more work 
should be done on this topic. 

In summary, we see several new research directions 
that can be opened by the present work. Firstly, sub¬ 
system DFT applications can surely benefit from the use 
of meta-GGA functionals which provide increased accu¬ 
racy with respect to GGAs at a lower computational cost 
than hybrid methods. In this context new meta-GGA XC 
functionals can be tested apart from TPSS. In second 
place, additional research can be done to develop more 
accurate semilocal approximations for the KED, which 
can be useful in the calculation of the non-additive XC 
energy. 

Finally, additional work can be foreseen to exploit the 
full power of the Laplacian level of theory, investigat¬ 
ing the use of Laplacian-level meta-GGA kinetic energy 
functionals in conjunction with similar approxima¬ 
tions employed in the non-additive XC term. In this way 
density and embedding errors in the kinetic and XC part 
can be expected to be more balanced. 
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Appendix A: Behavior of (20/9)r^^q in 

the tail of the density 


Let consider a spherical atom, where the density de¬ 
cays exponentially as p = Ae~°‘^, with r being the radial 
distance from the nucleus. In this case, 


I. 


(Al) 


q ar — 2 r- 

A similar expression can also be obtained for a Gaussian 
decaying density n = 

2ar^ 

q 


1 . 


2Q;r^ — 3 r-s-oo 

Thus, in the tail of the density q~ s^, and so 

20 

¥ 


20 ^ 
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FIG. 3: Plot of z = /t for the Ne atom and a non rela¬ 
tivistic atom with 1338 electrons, as computed with the exact 
T, T = , and r = r^. 


where we considered that the revAPBEK enhancement 
factor asymptotically behaves as a constant. Finally we 
find 


5/3 _ 3 
]—>oo 20/9 4 


(A4) 


Note, however, that this behavior is valid only at large 
distances that are not important in practical calculations, 
see Fig. |3l In valence and close tail regions z is instead 
very close to the exact one, as shown in Fig. [1] and Fig. 
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